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ABSTRACT 

A technique is presented for solving the inverse dynamics of flexible 
planar multibody systems. This technique yields the non-causal joint efforts 
(inverse dynamics) as well as the internal states (inverse kinematics) that pro- 
duce a prescribed nominal trajectory of the end effector. A non-recursive glo- 
bal Lagrangian approach is used in formulating the equations of motion as well 
as in solving the inverse dynamics equations. Contrary to the recursive method 
previously presented, the proposed method solves the inverse problem in a sys- 
tematic and direct manner for both open-chain as well as closed-chain 
configurations. Numerical simulation shows that the proposed procedure pro- 
vides an excellent tracking of the desired end effector trajectory. 


1. Introduction 

Accurate positioning and vibration minimization of flexible multibody systems have gen- 
erated considerable interest from the computational dynamics and controls communities. The 
advent of the new generation of very fast, lightweight robots and flexible articulated space 
structures has made the control of structural vibrations an important practical problem in the 
manufacturing and space industries. 

There is a large body of literature dealing with the forward dynamic analysis of flexible 
multibody systems, i.e., determining the resulting motion when the joint forces and external 
forces are given. Numerous approaches have been proposed that are either based on the mov- 
ing frame method or the inertial frame counterpart (see reference [1] and references therein.) 
Similarly, numerous control approaches have also been proposed for position and vibration 
control of flexible articulated structures (see reference [2] and references therein.) 

Solutions to the non-collocated control of flexible articulated structures have been 
presented in [3-6]. The so-called inverse dynamics joint actuation are non-causal or time- 
delayed joint torques (applied in negative time and future time) that are capable of positioning 
the end effector according to a desired trajectory. The importance of using the inverse dynam- 
ics approach to vibration control has been demonstrated recently in reference [7] where passive 
feedback and feedforward of the inverse dynamics torque were used to achieve an exponen- 
tially stable tracking control law that yields excellent end-point tracking of flexible articulated 
structures. In this paper, present a global Lagrangian approach to the solution of vibration 
minimization and end-point trajectory tracking. 
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2. Mathematical Formulation 

In order to describe the dynamic modeling let us consider a generic flexible body (Fig. 1) 
representing a component of a flexible articulated structure. The configuration of the multi- 
body system can be described by two sets of coordinates: the first set corresponds to the rigid 
body coordinates representing the location and orientation of the body axes, with respect to the 
inertial frame; and the second set corresponds to the so-called deformation coordinates or nodal 
deformations representing the deformation of the body with respect to the body axes. Using 
the aforementioned choice of coordinates, the location of an arbitrary point P in a planar 
deformable body / is given by 

i-' = R 1 ' + A‘ u* 0) 


where R‘ is the location of the origin of the body axes with respect to the inertial frame, u‘ is 
the location of point P with respect to the body axes, and A* is the rotation transformation 
matrix from the body axes to the inertial frame. In the three-dimensional case, the rotation 
transformation matrix is given by 


A‘ = 


2(0o + e, 2 ) - 1 
2(0,02 + 0 O 0 3 ) 
2 ( 0,03 - 0 O 0 2 ) 


2(0102 “ 0O®3) 
2(0o + 0 2 2 ) ~ 1 

2(0203 + 0O®l) 


2(0i Oj + 0002) 

2(0203 “ 000 1 ) 
2(0q + 0 3 2 ) - 1 


( 2 ) 


where the orientation coordinates are represented by four Euler parameters 0 q, 0,\ 0 2 , and 0 3 


which satisfy the following identity: 

£ (ei) 2 = 1 

* =0 


The vector u' can be decomposed into 

u‘ = u‘ + uj (3) 

where u‘ is the position vector of point P in the undeformed state with respect to the body 
axes, and uj is the deformation vector of point P with respect to the body axes. Differentiating 
Eq. (1) with respect to time yields the velocity vector of point P 

r‘ = R' + A' u* + A‘ u/ (4) 

where ( ) represents differentiation with respect to time. To separate the generalized coordi- 
nates, the second term on the right hand side of Eq. (4) may be written as 

A‘ u‘ = -2 A' u‘ E‘ 0‘ (5) 


where E‘ is a matrix that depends linearly on the Euler parameters and is given by 


E‘ = 


and u' is a 3 x 


— 01 00 03 — 02 

—02 — 03 00 ®1 

—0 3 02 —0, 0o 

3 skew-symmetric matrix given by 


u 


0 ~ U z M > 

u, 0 -u, 

-u y 0 


( 6 ) 


(7) 
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where u x , Uy, and u, are the coordinates of the generic point P with respect to the body axes, 
in the deformed configuration. 

The deformation vector uj can be expressed in terms of the nodal deformations by using 
a finite element discretization scheme 

u / = N 1 ' q/ (8) 

where N* is the shape function matrix and q| is the nodal deformation vector. Since the shape 
function matrix is time-invariant, the time derivative of the deformation vector becomes 

u; = N‘ q; (9) 

where q } is the time derivative of the vector of nodal deformations. Substituting Eqs. (5) and 
(9) into Eq. (4), we obtain the following expression for the velocity vector in terms of the rigid 
body coordinates and nodal deformation coordinates: 

r* = R‘ - 2 A 1 u‘ E‘ 0' + A' N' q } (10) 


Using Eq. (10) to describe the velocity vector of an arbitrary point P , the kinetic energy 
of body i can be expressed in the following quadratic form in velocities 


(ID 
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where the constant submatrices % and vn,r represent the total mass of the body and the con- 
sistent finite element mass matrix, respectively. The submatrix m 00 represents the moment of 
inertia of the deformable body about the origin of the body axes, and the submatrix 
represents the first moment of mass of the deformable body about the body axes. The subma- 
trices m f R and m/ 0 represent the kinematic coupling between the rigid body coordinates and 
the nodal deformation coordinates. 


The potential energy due to linear elastic strains in the material can be expressed in the 
following quadratic form in rigid body coordinates and nodal deformation coordinates 

( 12 ) 

where k is the conventional finite element stiffness matrix. The singularity of the stiffness 
matrix can be eliminated by imposing appropriate boundary conditions or by choosing vibra- 
tion modes that are consistent with the boundary conditions. 


P.E } = 
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2.1. Equations of Motion 


In order to unify the equations formulated in rigid body dynamics and structural dynam- 
ics, we make use of generalized coordinates which include rigid body coordinates and deforma- 
tion coordinates, hence 


I 



( 13 ) 
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where q' is the vector of generalized coordinates for body i, q‘ and q/ are the rigid bodv 

S f “d C r inmS ° f ^ '■ ^ eneiwor 


K.E . ' = — 
2 


q‘ T M‘ q' 


( 14 ) 


where M' is the overall mass matrix of body i. Similarly, the potential energy of the bodv due 
to linear elastic deformation can be expressed by ^ 


PE. 1 = — q ,; 
2 M 


K' q' 


( 15 ) 


where K 1 is the overall stiffness matrix of body i. 

When reference coordinates such as those described in this paper are employed in multi- 
i y syst ^ f s ’ generalized coordinates are not independent because the motion of specific 
points m different bodies are related according to the type of mechanical joint that connects the 

flCXib " mecha " icaI s >™ ras - <l=f»rmSion of a c^nem 
affects the configuration of adjacent components. The inteniependence of the generalized mor 
dmates are expressed by a vector of kinematic constraint equations, such as § “ 

<D(q,r) = 0 

. - u (16) 

wJri q - ,s . the l ° tal v u cc , tor of system generalized coordinates, t is time, and O is the vector of 
^ '" dc Pf nden ‘ holonomic constraint equations. For example, if a revolme iota co^ecm 


[r‘ + A‘ u>] - + A' U^J = 0 


( 17 ) 

In generai, holonomic constraints can also be explicit functions of time as well as generalized 
traS^T’ ^ m CaSC ° f imP ° Sing thC C °° rdinates 0f me end-effector to follow a dtol 

lions o/mLonSSeT to ^ * ***** ^ COnstrained eoordinates, the system equa- 


M(q) q + C q + K q + <D q T X = Q t + Q^q) 


( 18 ) 

rcsrectitdy Vtoe vermin fT StCm m3SS ’ , Sy ? em dam P in g 30(1 system stiffness matrices, 
respectively, X is the vector of Lagrange multipliers associated with the constraints <D is the 

constraint Jacobian matnx, Q e is the vector of applied external forces, and Q is the quadratic 

velocity vector. The quadratic velocity vector contains the centrifugal forces and Coriolis 

geSS SiM«.‘ h£ difreren,iafi0n of ki « ic expression with mspec, ,o ihe 

quations (DAE s) that have to be solved simultaneously. The solution to the inverse dvnam 
ics problem requires a forward dynamic analysis within an iteration process We solve the for 
ward dynamrcs problem by using the augmented Lagrangian penal, yTtoatl m The au« 
agrangian penalty formulation obviates the need to solve a mixed set of differential- 
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stSntTVpp^g ^au^ented teams’ an* ° f equations t0 account for the con- 

result in the following equation: ^ gI fom, uIation to Eqs. (16) and (18) will 
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( 20 ) 

22. Inverse Kinematics and Inverse Dynamics 
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Eq. (23) results in 

m ff + c // <1/ + k // q / = Gy m e/ qy + F,(X,q r ,q r ,q r ,q f ,q f ) (25) 

where Fi is a force vector that includes the inertial terms, reaction terms between contiguous 
bodies, and quadratic velocity terms. 

The inertial coupling submatrix m 0 y can be decomposed into a sum of a time-invariant 
matrix and a time-varying matrix 

m 0/ = m e/ + m e/ (26) 

where m 0 y and m 0 y arc the time-invariant part and lime-varying part of m 0 y , respectively. 
This decomposition is essential for the iteration process needed to obtain the solution as 
explained below. Substituting Eq. (26) into Eq. (25), we obtain the inverse kinematics equa- 
tion of motion for the nodal displacements q y : 

m ff q / + c // q/ + k // q / = F(X,q r ,q r ,q r ,q f ,q f ,q f ) (27) 

where 


m ff 


= m ff - Gy m£y 


(28) 


The problem statement for the inverse kinematics is that of finding the internal states q , 
so that the end-point coordinates characterized by a subset of the rigid body coordinates q r fol- 
low a prescribed trajectory. The mass matrix myy is nonsymmetric and it is precisely the non- 
symmetry of the mass matrix that produces internal states which are non-causal with respect to 
the end-point motion. Eq. (27) is a nonlinear differential equation in the variable qy, As 
explained below, Eq. (27) is solved iteratively in the frequency domain to yield the nodal 
deformation vector qy that is non-causal with respect to the end-point motion. 

In the frequency domain, Eq. (27) can be written as a set of complex equations for a par- 
ticular frequency GJ 


m " + m ' F k// 


qy(ffi) = fxw) 


(29) 


where q^GJ) is the Fourier transform, of q 7 (r) and P(GJ) is the Fourier transform of F(r). Eq. 
(29) is based on the assumption that tjy(r) and F(r) are Fourier transformable. This assumption 
is valid for slewing motions which are from rest to rest. The nodal acceleration vector (^(55) 
can be obtained directly from Eq. (29) for each frequency G5. The leading matrix of Eq. (29) is 
a complex regular matrix that is invertible for all frequencies except for 55 = 0. However, for 
S5 = 0, the system undergoes a rigid body motion determined only by the invertible mass 
matrix rriyy . The nodal accelerations in the time domain may be obtained through the applica- 
tion of the inverse Fourier transform, i.e., 


q'/(0 = ^ J ^(C5)e ,cy dm 

— oo 


(30) 


Once the non-causal nodal accelerations are known, Eq. (22) can be used to explicitly 
compute the non-causal inverse dynamics joint efforts that will move the end effector accord- 
ing to a desired trajectory. We note, however, that the inverse dynamics torque and internal 
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states given by Eqs. (22) and (27), respectively, depend on the Lagrange multipliers and rigid 
body coordinates, which in turn depend on the internal states and the applied torque. Moreover, 
the rigid body coordinates and Lagrange multipliers are different from their nominal values 
when the components of the multibody system are flexible. Therefore, a forward dynamic 
analysis is required to obtain an improved estimate of the generalized coordinates and Lagrange 
multipliers given the torques computed previously using nominal values of rigid body coordi- 
nates and Lagrange multipliers. In order to ensure that the iteration process converges to obtain 
the joint efforts that will cause the end-effector to follow the desired trajectory, the forward 
dynamics analysis is carried out with the additional constraint that the coordinates of the end- 
point follow the desired trajectory. These additional constraints have corresponding Lagrange 
multipliers which act as correcting terms to the joint efforts that have been previously calcu- 
lated. 

To summarize, the procedure for obtaining the inverse dynamics solution for flexible mul- 
tibody systems involve the following steps: 

Algorithm: 

1. Perform a rigid body inverse dynamic analysis to obtain the nominal 
values of the rigid body coordinates q r and Lagrange multipliers A.. 

2. Solve the inverse kinematics equation represented by Eq. (27) 
to obtain the time-delayed nodal accelerations . 

3. Compute the inverse dynamics joint efforts Q i0 using Eq. (22). 

4. Perform a forward dynamic analysis using Eq. (19) to obtain new 
values for the generalized coordinates and Lagrange multipliers. 

5. Repeat steps 2 through 4 until convergence in the inverse dynamics 
torques is achieved. 


3. Simulation Results 

We present in this section the results of numerical simulations that verify the procedure 
discussed above. First, we apply the global Lagrangian approach to an open-chain flexible mul- 
tibody system and compare the results with those obtained by the recursive Newton-Euler 
method [5] to test the validity of the proposed procedure. Next, we present the results of the 
application of the global Lagrangian approach to a closed-chain flexible multibody system to 
determine the inverse dynamics torque that will produce the desired motion at the end effector. 

3.1. Open-Chain Multibody System 

Fig. 3 shows a two-link flexible multibody system in the horizontal plane. The end-point 
of the second link is specified to move along the x-axis according to the acceleration profile 
described by Fig. 5, which corresponds to an end-point displacement of 0.483 meters along the 
x-axis. The geometric and material properties of the links are: 
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First Link: 

Length: 0.66 m 

Cross-section area: 1.2 x 10 -4 m 2 

Cross-section moment of inertia: 2 3 x 10 -10 m 4 
Second Link: 

Length: 0.66 m 

Cross-section area: 4.0 x 10 ~ 5 m 2 
Cross-section moment of inertia: 8.5 x 10~ 12 m 4 

The two links share the following properties: 

Young’s modulus: 14 GPa 
Mass density: 2715 kg/m 3 


gian method is superim^sed oTthe 016 gIobal La g ran - 

sive Newton-Euler method The inverse rivnimLc ^ ue e determined by the recur- 
pu^ by me ,wo aforcmentimied SMTSMU'E Z * B 'JTf? ^ 

“ me »me ^ “ 

3.2, Closed-Chain Multibody System 

two jotatthicTtre t°^aS,man!’ltr' 11 'r ly sys T made " p of four »««>"= Unis wim 
■be multibody system Is ^n “°" T C '° **, 8ro,,nd - As in open-chain case, 

neglected. The dcsSmd tmjec oTof oi S i" a <?* vi » eff “ ™ 

and y axes. The x and yiommneni L accel2n„ ^ Wllh rcspect 10 ““ * 

acceleration profile shown in Fig. 8 The four links share ihVcut 5 315 speclfle<1 10 follow the 
properties: ■ ■ re tbe following geometric and material 


Length: 0.60 m 

Cross-section area: 4.0 x 10 ~ 5 m 2 
Cross-section moment of inertia: 8.5 x 10 -12 m 4 
Young’s modulus: 14 GPa 
Mass density: 2715 kg/m 3 
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body l„r q „e produces res, dual an^'S 
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residual angular rotations. As a matter of fact, it has been observed in the simulations that the 
rigid body torques produced residual vibration in all the nodal deformations while the inverse 
dynamics torques eliminated the residual oscillation. Furthermore, the inverse dynamics torques 
produced nodal deformations which exhibit non-causal characteristics with respect to the end- 
point motion. Fig. 12 shows a comparison of the vertical tip error at joint 5 obtained by a 
feedforward of the inverse dynamics torque with the tip error resulting from a feedforward of 
the rigid body torque. This figure shows that the inverse dynamics torque provides an excel- 
lent tracking of the desired end effector trajectory whereas the rigid body torque again induces 
substantial vibration on the end-point motion. 

4. Conclusion 

A global Lagrangian approach for the inverse dynamics of flexible multibody systems has 
been presented. The procedure is capable of solving for the inverse dynamics torque profiles of 
both open-chain and closed-chain flexible multibody systems in a unified and systematic 
manner. The method is found to produce an excellent tracking of the desired trajectory of the 
end effector. In a future paper, we will address the inverse dynamics problem for flexible mul- 
tibody systems undergoing motion in three dimensions. New problems arise in the three- 
dimensional case, since the actuating torque vectors have directions which are time-varying and 
nonlinear functions of the rigid body coordinates, as contrasted with the planar case where the 
applied torque vectors have directions fixed perpendicular to the plane of the multibody system. 
In addition, to be able to perform the inverse kinematics and inverse dynamics analyses, addi- 
tional actuation at the joints may be necessary. 
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